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ABSTRACT 

We derive expressions for the tidal field exerted by a spherically symmetric galaxy having an 
extended mass distribution, and use our analysis to calculate tidal heating of stars in a globular 
cluster or a satellite galaxy. We integrate the tidal force over accurate orbits in various models for 
the primary by using the impulse approximation, coupled with adiabatic corrections that allow for 
time variation of the perturbation. We verify the results with direct N-body simulations, and also 
compare them with the conventional straight path approximation. Heating on highly eccentric 
orbits dominates, as adiabatic corrections strictly prohibit energy changes on low eccentricity 
orbits. The results are illustrated for globular cluster NGC 6712. For the orbital eccentricities 
higher than 0.9, the future lifetime of NGC 6712 is less than 10^" yr. Our analysis can be used 
to study tidal interactions of star clusters and galaxies in a semi-analytical manner. 



Subject headings: globular clusters: general - galaxies: interactions 



1. Introduction 

Tidal forces are very important in stellar systems. They determine the dynamics of interactions between 
galaxies in clusters, between a dwarf satellite and a primary galaxy, and between a star cluster and its host 
galaxy. The most dramatic tidal perturbations occur in fast and close encounters. When the duration of the 
encounter is shorter than the characteristic dynamical time of the satellite, such an interaction is referred to 



as a tidal shock (Spitzer 1987) 



Many previous semi-analytic and numerical studies of tidal shocks considered the primary galaxy as 
a point-mass ( Spitzer 1958 ; Richstonc 1975 ; Knobloch 1976 ). Such an approximation is violated when a 
satellite is well within the limits of the galaxy or a globular cluster passes near the galactic nucleus. More 



recent work included the extended structure of the perturber (King and de Vaucouleurs models in Aguilar 
fe White] 1985, 1986; isothermal sphere potential in ph fc Lin 1992| ; ph, Lin fc Aarseth 1992| , 1995; and 



synthetic models in Johnston et al. 1995; Kroupa 1997), but most results were obtained through numerical 
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simulations. Weinberg (1997) and Murali & Weinberg (1997a-c) used the linear perturbation theory and 



semi-analytical Fokker-Planck calculations in their study of globular clusters and elliptical galaxies. The 
latter work is complex and cannot be easily used for quick estimates. In this paper, we obtain expressions 
for the tidal field of an arbitrary spherically-symmetric system and apply them to estimate the amount of 
heating of the satellite stars produced in the interaction. The resulting expressions are simple and ready-to- 
use. 

The paper is organized so that we move from the simplest and least accurate treatment to the more 
complex and more accurate one. We present analytical calculations of the tidal field and energy changes for 
straight-path orbits in §^ and for eccentric orbits in §|^. First, we use the impulse approximation and then 
relax it by allowing for the motion of stars during the perturbation. Some tedious integrals are relegated to 
the Appendix. Then, in we apply our theory to calculate the heating of globular cluster NGC 6712 in 
the Galaxy. We discuss and compare our method with previous work in ^ In we summarize our results 
for semi-analytic studies of tidal interactions. 



2. Tidal perturbation in straight-path encounters 

Consider an inertial reference frame with the origin at the Galactic center. Let Ro(t) be the trajectory 
of the center of the satellite, and r the distance of a member star from the satellite center. Then the 
radius- vector of a star is R = Ro(0 + ^- The equations of motion are 

R- (U 

or 

+ ^ - (1^) - (^) , - (^01;;) ^ - - 

where and <I>g are the potentials of the satellite and the galaxy, respectively. We will drop the subscripts 
henceforth. The trajectory of the center of mass of the satellite is determined by the leading term, Rq = 
—d^c/dRo, so, in the tidal approximation, we have 

^ dr dRodRo ' 

The potential of a spherically-symmetric system is 

$G = _ r p{r')r'dr', (4) 

R Jr 

which produces a tidal force per unit mass 

d^^G_\ _ GM q 
dRadRo) 

where AIq is the total mass of the galaxy, and /^(i?) is the normalized mass profile 



Ftid = - ( an ar. I = ^53" [(^M " M)(n ' r)n - pir] , (5) 



and /i is defined by 
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The direction to the center of the satelHte is denoted by n = Ro/i?o- 

The treatment so far in this section is general. We will now temporarily restrict ourselves to the impulse 
approximation and straight line orbits for simplicity in the computation. In the impulse approximation, we 
neglect internal motions of the satellite stars during the encounter. The corresponding change in the star's 
velocity is the integral 

Av = / Fud dt (8) 



over the duration of the interaction. This time-varying perturbation leads to an increase of the random 
motion (or heating) of stars in the satellite. When the phase space is well mixed, the averaged energy change 
of stars with the initial energy E is quadratic in perturbation: 

{AE)E = {\{Avf). (9) 

We have also neglected the self-consistent reaction of the satellite potential as it settles into a new post-shock 
equilibrium. Effects of the potential readjustment are studied elsewhere ( Gnedin fc Ostriker 1998). 

To calculate the velocity kick, Av, we need to specify satellite's orbit. Since the tidal force is short-range, 
only the part of the orbit within a few distances of the closest approach (perigalacticon Rp) is important. 
Therefore, the orbit is conventionally assumed to be a straight path, Ro(i) — {Rp,Vpt,Q)^ where Vp is the 
orbital velocity at perigalacticon. In §^ we consider eccentric orbits of the satellite. 

Integrating equation (^) from t = — oo to oo, we find 

= -BtTT {(3^0 -Ji- h)x, (2/o -h- 3Jo + Ji)y, -hz} , (10) 



RlVp 



where 



/o(i?p) EE KRpC) ^2 _ -^^1/2 , (11a) 

h{Rp) ^ KRvC) ^2 (^2 _ 1)1/2 ^ (lib) 

f°° dC 

MRp) ^ A^(^pC) ^i^^y^YyiTJ, (11c) 

r°° dc 

MRp) EE / fiiRpC) ^4.^2 .M/2 - (lid) 



Notice also that 

■^'^^"dltp- ^^^^ 

Averaging over an ensemble of stars in a spherically-symmetric satellite, we have (x^) ~ {y^) = {z^) = 
The heating term per unit mass is therefore 

Xst = i [(3Jo- Ji-/o)' + (2/o-Ji-3Jo + Ji)2+/o']. (14) 

These results generalize the usual expressions for a point-mass (e.g. Binney fc Tremain^ 1987) to the case 
of an arbitrary spherical density profile. In particular, the function XstiRp) is the correction due to the 



-4- 



extended mass distribution of the primary galaxy. If all galactic mass was concentrated at its center, then 
yu — 1 and il — everywhere, and we have 



/n = 



= 1, 



Jn = 



(15) 



and II = Ji = 0. The factor Xst becomes unity in this case, as expected. The subscript "st" reminds us that 
this factor is valid only for the straight path approximation. 



The second order energy diffusion term to the lowest order in perturbation is 



{AE^)e = ((v • Av)2) = 



2GMo 
R^Vp 



2r2 V^l + Xr,v) 



(16) 



where Xr,v is the position-velocity correlation function, which takes values from -0.25 to -0.57 (see Gnedin 
fe Ostriker 1998). This can be simply rewritten as 



(AE). 



(17) 



Higher order terms, {AE"), are negligible for tidal shocks ( ^pitzcr 1987 ). For example, both n = 3 and 
n — A terms are quadratic in the perturbation energy, AE / E <C 1. Therefore, only the first two terms are 
important in the Fokker-Planck equation (see §0) . 



2.1. The Hernquist model 



For illustration, we consider several examples of mass distributions, the Hernquist model (Hernquist 



1990| ) and power-law density profiles. The Hernquist model closely resembles the i?^/^ law for spherical 



galaxies and gives simple analytical expressions for the potential-density pair: 

$ = -^ p= (18) 

i? + a' ^ 27ri?(i? + a)3 ' ^ ' 

where Afo is the total mass of the galaxy, and a is the scale length. The normalized mass profile is 



The calculation of the integrals ( pi] ) is given in Appendix A. The results (eq. [A9|) are plotted in Figure |l| 
versus the dimensionless parameter a = a/ Rp. 

For large impact parameters {a <C 1), both Jq and Jo approach constant values {Iq = 1, Jo = 2/3), 
while II and Ji decrease linearly with a. As a result, the satellite is compressed in the vertical direction, 
stretched in the direction of the perigalacticon by the same amount, and unaffected in the perpendicular 
direction. All this is expected in the point-mass approximation. At very small impact parameters (a 3> 1), 
when the satellite penetrates within a scale length of the center, all correction factors decline approximately 
as a~^. This asymptotic behavior is due to the mass of the model increasing quadratically with the radius 
at small r. 

In Figure ^, we plot the absolute value of the components of the velocity change for the Hernquist 
model: 

Xx = 3 Jo - Ji - Jo (20a) 
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Xy = 2/0 - Jl - 3 Jo + Jl (20b) 
Xz = h (20c) 

Note that Xx and Xz are the velocity changes in the x and z directions relative to the point-mass ap- 
proximation. The velocity change in the y direction is negligible for any a; its typical value is Xy ^ 10^"*. 
In the point-mass approximation, Auj, = 0. Also, we plot for comparison a naive approximation, which 
assumes the galactic mass inside the perigalacticon Rp to be a point-mass. The correction factor in this case 
is XxQ — X20 — Xq ^ l^iRp)- It underestimates the vertical compression, but overestimates the effect in the 
radial direction. 



2.2. The power-law models 

Next, we consider a general power-law distribution truncated at a limiting radius Rmax, 

p{R)^{ R<Rrnax (« < 3) ^^1) 

t 0, R > Rm ax 

Here fc is a constant which determines the total mass of the model, Afo = ^z^^^"- We assume that the 
perturbed satellite is well within the extent of the galaxy, so that truncation of the model does not affect 
our results. The normalized mass profile is 

' m = iS-n)fi{R). (22) 
The evaluation of the tidal integrals (eqs. 0|) is straightforward: 

Note that the integral Jo converges only for n > 1. Using the identity T{x + I) = xr{x), we find n Jq = 
(n-I)Jo. 

The change of the stellar velocity is 

where the y-component vanishes because of the symmetry of the trajectory. The correction to the energy 
change is 



Xst 



Rma:i 



r(f) 



We remind the reader that this result holds for the power- law indices 1 < n < 3. 

Notice the change in sign of Ati^; as the index n passes through the value n = 2. For the special case 
of an isothermal sphere (n = 2), the x-component of the velocity change vanishes identically, and the only 
resulting change is compression in the vertical direction: 

GMo TTRp TT^ Rl 

^^^~mv~R^^ X-t(-Rp) - (26) 
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3. Tidal perturbation on eccentric orbits 

The straight-path orbit is a convenient approximation, but in real galaxies, satellites move on eccentric 
orbits. Effects of tidal shocks are most prominent in the vicinity of the perigalacticon. The results obtained 
in the previous section should nevertheless give us insight into the actual heating of satellites on eccentric 
orbits. In this section, we determine the true orbits in our examples of galactic models and investigate how 
tidal heating depends on parameters of the orbit. 

In a spherically-symmetric potential, the orbit of a satellite is confined to a plane and can be charac- 
terized by two independent parameters. Let us choose the perigalactic distance, Rp, and eccentricity, e. as 
such parameters. They in turn determine the energy and angular momentum of the orbit. We employ polar 
coordinates, in which the position angle 9{t) can be used as a time variable. This angle is defined such that 
6 = Q aX R = Rp. In the coordinate frame where the orbit is in the XY plane, the position vector of the 
satellite is n(t) = xcos^(f) -l-ysin^(t). 

Using the impulse approximation, we integrate the tidal force along the orbit to obtain the expected 
change in the stellar velocity: 



.T/2 

Av= / Ftidrft, (27) 

J-T/2 

where T is the orbital period. 



Changing the variable of integration from time, t, to angle, 9, using the identity dt = {R^ / J) dO, where 
J is the orbital angular momentum, we obtain 

Av = dO ^ [(3/. - m) (n • r)n - ^ , (28) 

where a is the parameter of the galactic model {Rmax in case of the power-law models), and j is the 
dimensionless angular momentum, defined by J = {GMq a^^'^j . The integration over the position angle 
extends to the maximum angle 9„i corresponding to the apogalactic point of the orbit. Thus, the tidal force 
is integrated starting from the apogalacticon at 6* = —9m, approaching the perigalacticon at = 0, and 
going to the next apogalacticon at 9 = 9m- The maximum value of the position angle ranges from 9m = 
in the Kcplerian potential to 6*™ = 7r/2 in the harmonic oscillator potential. A very accurate procedure for 
calculating the orbits is described in Appendix B. 



Rewriting the velocity change into components, we have 

V a3 y j{a,e) 



f CM \ ^^"^ 1 

Av= ^ {{Bi-B^)x,{B2-B:i)y,-B3z}, (29) 



where 



B, ^ I ' ^J^^^0^cos^9d9 (30a) 



B2 

Bs 



' 3Mi^)-M(i^) 
R/a 

= r ^d9. (30c) 
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The cross-terms (cx xy) vanish because of the symmetry of the orbit, R{0) = R{—9). The energy changes 
of stars with the initial energy E are 



6f{a,e) 



We have now added adiabatic corrections, Ai{xt) and A2{xt), which allow for the motion of stars during 
the perturbation. Adiabatic corrections reduce the amount of heating as the result of the conservation 



of stellar adiabatic invariants. Gncdin & Ostriker (1998) provide a detailed discussion of the adiabatic 
corrections in the case of disk shocking. The main results can be summarized as follows. The change of 
the stellar energy inferred from the N-body modeling of tidal shocks is described to a good accuracy by the 
product of the impulsive value and the adiabatic correction of the form 

A,(xO = (l + a.?)-^% (33) 

where Xt = lot. Here ijj{r) is the orbital frequency of stars in the satellite, and t is the effective duration of 
the shock. The exponents 7^ depend on the shock duration relative to the half-mass dynamical time of the 
satellite, tdym and vary from (71 = 2.5, 72 = 3) for r <^ tdyn to (71 = 72 = 1.5) for r ^ 4:tdyn- The adiabatic 
correction becomes increasingly small in the satellite core, conserving stellar actions and energy. 

The simulations were done for the King model of the satellite, with the concentration parameter c = 0.84. 
We have checked that expression ( |3^ ) remains valid for most parts of the more concentrated model with 
c = 1.5. However, our adiabatic corrections slightly underestimate the energy change in the core. It is 
possible that equation (BSf) breaks for very concentrated satellites and thus should be applied with caution. 



See also Weinberg (1994a-c) for the direct application of the perturbation theory formalism. 



The adiabatic correction was calculated under the assumption that the perturbing force varies with 

.2 / 2 . 

time as a Gaussian function, i^tid oc e~ '"^ . This assumption does not strongly limit the applicability of the 
results since numerical simulations show that significant energy changes occur only close to the perigalacticon, 
where almost any tidal force has similar behavior. As an example, consider a Keplerian orbit of eccentricity 
e. Near the perigalacticon, the angle 9 <ti 1 and the vertical component of the tidal force varies in agreement 



with the Gaussian: Fzz ~ ^"22,0(1 ~ ci9^) « i^2z,o(l ~ C2t^), where ci and C2 are some constants. See §4.1 
for more discussion of the temporal structure of the tidal force. 

The adiabatic correction depends not only on the stellar frequency, uj, but also on the orbital eccentricity 
of the satellite, through the parameter r. On low eccentricity orbits, the effective duration of the "shock" 
becomes so long that adiabatic corrections prohibit the changes of the stellar energy. On circular orbits, the 
effects of tidal shocks are greatly reduced. We will return to adiabatic corrections in when we consider 
the evolution of globular cluster NGC 6712. 



3.1. The Hernquist model 

It is instructive to compare the results for the satellites on eccentric orbits with those on the straight path. 
For the Hernquist model, we perform the integration of equations ( |30| ) numerically for a grid of perigalactica 
and eccentricities. We now rewrite equation (pT|) for (Ai?) similarly to equation (O), extracting the factor 
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Xecc which is now a function of the two parameters, a = a/ Rp and e. For clarity, the factor Xecc{ce, e) does 
not include the adiabatic correction Ai(xt). 

Figure m shows the ratio Xecc/Xst for the orbits with various eccentricities. For highly eccentric orbits, 
e < 1, the ratio varies shghtly with a, from Xecc/Xst ~ 1-8 for a <C 1 to Xecc/Xst = 1 for a 1. The latter 
limit corresponds to the most eccentric Keplerian orbit where both approaches must agree. However, as the 
eccentricity decreases, the ratio Xecc/Xst varies with a much more strongly. While at a ^ 1, the ratio is 
highest for e = 0, the opposite is true for a 1. 



3.2. The power-law models 

The power-law models are scale-free and allow further simplification of equations (|3^). We have 
Bi{a, e) = q;"~^ Ci{e), p{a, e) = a"^** /(e), and therefore 

C AT 

(AE) = —J^r'a"xecc{e)A,{xt), (34) 

Xecc - . (35 j 

The remaining integrals depend only on the eccentricity e. After some algebra, we obtain 

Ci = n / ( — ) cos'^ ed9, (36a) 



C2 ^ n I I ^ ) sin^ede, (36b) 



C3 = / ( — ) de. (36c) 



Note that in the harmonic oscillator potential (n = 0), the first two integrals are identically zero: 
Ci = C2 = 0. The third one can be done using simple analytic orbits allowed in this potential (eq. [B2|): 

1 + e TT^ 

C3 = TT , Xecc = (for = 0)- (37) 

1 — e 2 

Surprisingly enough, the change of the stellar energy depends neither on the location in the galaxy, a, nor 
on the orbital eccentricity, e. The satellite is equally compressed in all directions: 

Av = I— — I n{-x,-y,-z}, (38) 

^'■max ^ 



(AiJ) = ^f.r'A,{x,) (forn^O). (39) 



In the isothermal sphere potential (n = 2), the integrals are also done straightforwardly: 

Ci = 2e',„ + sin 26l„„ (40a) 

C2 - 20™-sin20„, (40b) 

C3 = 29^, (torn = 2). (40c) 
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For the n ~ 1 model, the integrals Ci{e) are done numerically. Figure ^ shows the factors Xecc for all 
three models. For comparison, we also plot Xst for the isothermal sphere model. Heating on eccentric orbits 
is larger than that on the straight path, except for e = 1, where the two are equal. However, the plot does 
not include adiabatic corrections, which would have effectively reduced heating on low-eccentricity orbits. 



4. A Case Study: NGC 6712 

We now illustrate the theory developed in by considering the evolution of a globular star cluster, 
NGC 6712. This cluster is presently at the distance of 3.8 kpc from the Galactic center, close enough to 
experience tidal shocking by the Galactic bulge. We use the observed parameters to construct cluster's orbit 
in the bulge potential and to integrate the expected energy change. Then, we perform N-body simulations 
to verify the result. Finally, we consider various orbital eccentricities allowed by the cluster orbital energy 
and illustrate the eccentricity effect on the resulting heating. 



4.1. Constructing the orbit 

Full three-dimensional kinematic data for NGC 6712 are not yet available. Instead, we use the results of 



Monte-Carlo simulations by Gnedin & Ostriker (1997) for the Ostriker-Caldwell (1983) model of the Galaxy 



with the isotropic velocity distribution of the globular cluster system. We draw the two unknown components 
of the orbital velocity consistent with the chosen kinematic model and integrate the orbits for 10^° yr. For 
NGC 6712, the median perigalactic distance is Rp — 0.73 kpc and the velocity at perigalacticon is = 396 
km s^^. The eccentricity of this orbit is e = 0.69. 

We approximate the Galactic bulge by a Hernquist model with the scale-length a = 0.6 kpc. The total 
mass of the model required to reproduce the observed velocity Vp is Mq = 3.3 x 10^" Mq, close to the 
observed mass of the bulge. 

Figure || shows the vertical component of the tidal force exerted on the cluster along a single orbit in 
the Galaxy. The force peaks at the perigalacticon. The dashed line shows that the Gaussian fit to the actual 
force, Ffit oc e~* ^"^ , is accurate for the most important part of the orbit. (A similar parameterization was 



employed by Johnston, Hernquist fc Weinberg 1998 in their numerical study of tidal heating.) The width of 
the Gaussian, or the effective duration of the shock, is r w 5.3trfj,„, where tdyn is the half-mass dynamical 
time of the cluster. This timescale is long enough that most stars in the cluster are protected by adiabatic 
invariants, even for Rp < 1 kpc. Using the calculated value of r, we estimate adiabatic corrections to the 
energy change (cf. eq. p3|) and verify them with N-body simulations described next. 



4.2. Self-Consistent N-body simulations 



We can now simulate the cluster as a system of N point-mass particles and run it along the chosen orbit. 
At the end, we calculate the changes of the stellar energies and compare them with our analytic theory. 



We use a Self-Consistent Field code ( [Hernquist fc Ostriker 1992| , [Cnedin fc Ostrikei] 1998) with N = 10^ 
particles. The code approximates the true gravitational potential and density of the cluster with a finite 
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series of basis functions. In this work, we take n = 10 radial functions and I = 6 spherical harmonics. 

The initial cluster realization is a King model with the concentration parameter c = 0.84. It is sufficiently 
close to the observed concentration of NGC 6712, c 0.9. The code employs units in which G = 1, M^i = 1, 
and i?c = 1, where M^i and Rc are the mass and the initial core radius of the cluster. The physical units 
are fixed by the observed values, = 2.5 x 10^ Mq and Rc — l-Q pc. 

In this simulation, the gravitational potential of the cluster is kept fixed. The self-consistent response 



of the system to tidal perturbations is studied elsewhere (Gnedin & Ostriker 1998) and is not important for 
our discussion. We run the cluster along the orbit, starting from the apogalacticon, where the tidal force is 
minimal (cf. Figure going to the perigalacticon, and then moving to the next apogalacticon. The whole 
simulation lasts for about 200 dynamical times, tdym of the cluster, with a short enough time step to ensure 
accurate integration of stellar orbits; Ai = 0.01 tdyn. We have run the simulation on the SGI Origin 2000 
supercomputer at the Princeton University Observatory. 

Figure ^ shows the first and second order energy changes at the end of the simulation. The particles are 
grouped in equal size bins, ranked by their initial energy, and the mean values, (Ai?) and (AE^), and their 
standard errors are calculated for each bin. The data are compared with our analytic result (eqs. [^ , ^ ), 
wherein adiabatic corrections (eq. p3| with 71 = 1.4, 72 = 1.5) are appropriate for the long duration of the 



shock, r, estimated in §4.1. Had we used the simplest theory, presented in §|^, for straight line orbits without 
adabatic corrections, we would have overestimated the total energy change significantly. By calculating the 
correct orbits and employing adiabatic corrections, we find a very good agreement with the simulations. 
This excellent agreement encourages one to use the semi-analytical theory (eqs. |^l[ - |3^ ) to study the effects 
of tidal shocks on the evolution of star clusters and satellite galaxies. 



4.3. Effects of orbital eccentricity 

Finally, we demonstrate how the eccentricity of the orbit affects the resulting heating. We keep the 
orbital energy of NGC 6712 fixed and consider all allowed eccentricities. Since these two parameters {Eorb 
and e) specify the orbit, the perigalactic distance Rp varies accordingly with e. Increasing eccentricity 
decreases Rp and vice versa. Also, low eccentricity orbits have longer shock durations, r, which leads to 
significant adiabatic corrections. 

Figure ^ shows the energy change at the half-mass radius of the cluster, calculated from equations 
(^,^^. The value of {/S.E)h falls dramatically with decreasing eccentricity, in part because Rp increases 
as e decreases and largely because of the adiabatic corrections. As expected, there is essentially no heating 
on nearly circular orbits. Also, we calculate the median energy change, {AE)h, for two more cases, when 
the orbital energy is a half and twice the observed value. The more bound orbits tend to have a shorter 
timescale r and correspondingly larger energy input than the less bound orbits. 



Discussion 



Our analysis applies to weak gravitational perturbations, tidal shocks. A significant evolution of a 
satellite or a star cluster occurs over a number of successive shocks. In this respect, the current work differs 
from that of Aguilar fc White (1985,1986) who studied the effects of strong interactions between two similar 
mass galaxies. Their work showed that the impulsive approximation is adequate for the most regions where 
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a significant mass loss and energy changes occur. 
I I 

Aguilar & White emphasized that particles escaping the satellite after the perturbation should be 
excluded from the analysis of the remaining structure. In the case of tidal shocks, few stars leave the 
satellite immediately. Instead, they slowly drift in the energy space (see Figure ^ and eventually escape 
through the tidal boundary. 

The results of our calculations are useful for semi-analytical studies of tidal interactions between large 
galaxies and their dwarf companions or star clusters. For example, Gncdin & Ostriker (1997) estimated the 



amount of heating of globular clusters due to the Galactic bulge by inserting equations (|ll|, |13; p^ ) into their 
Fokker-Planck code. Assuming the spherical symmetry of the cluster, the code solves the Fokker-Planck 
equation for the distribution function of stars in energy space. Since the amplitude of perturbations is small. 



the first two terms of the expansion of the Master equation (for example, [van Kampcn 198l| ; [Gardiner 1985| ) 
are both linear in the perturbation energy and lead to the Fokker-Planck formulation. The higher order 
terms are negligible and can be ignored. This approach is much faster than a full N-body simulation and 
allows a direct comparison with analytical models of star clusters. A new detailed study of the globular 



cluster evolution including tidal shocks is presented in Gnedin, Lee, & Ostriker (1998) 



6. Summary 

We have calculated the tidal field of a spherical mass distribution. Using the impulse approximation, 
coupled with adiabatic corrections, we estimated the amount of heating of stars in the satellite during a tidal 
shock. We considered several mass distributions for the host galaxy, a Hernquist model and three power- 
law profiles, and explored various satellite orbits. Although we restricted ourselves to spherically-symmetric 
galaxies to simplify the analysis, our results should also hold for non-spherical systems with a modest amount 
of flattening. 

For the Hernquist model, we found that heating is enhanced by a factor of few on eccentric orbits 
compared to the conventional straight path approximation. However, adiabatic corrections prevent any 
heating on low-eccentricity orbits. In the harmonic oscillator potential, the satellite is equally compressed 
in all directions. In the isothermal sphere potential, heating on eccentric orbits is also larger than that in 
the straight path approximation. 



To illustrate the analytic results (eqs. |31-^), we considered the example of globular cluster NGC 6712. 
For a fixed orbital energy of the cluster, the heating is much larger on highly eccentric orbits than on low- 
eccentricity orbits because of the adiabatic corrections. N-body simulations confirm the validity of our 
method, which becomes a powerful tool in calculating the tidal heating of satellites of larger galaxies. This 
semi-analytic method can also be used for studying the evolution of giant star clusters in galactic nuclei or 
the disruption of dwarf satellites in groups of galaxies. 



We would like to thank Simon White and the anonymous referee for useful comments. This work was 
supported in part by the NSF under grants AST 94-24416 and ASC 93-18185, and by the Presidential Faculty 
Fellows Program. 
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A. Calculation of the tidal integrals for the Hernquist model 



We evaluate here the tidal integrals (eq. jll]) for the Hernquist model. The normalized mass profile is 
given by equation (p^). Denoting a = a/Rp, the first integral is 



loia) 



dc 



1 ic+a)^ cHe~iy^'' 



Substituting x = ((^ + a) ^, we have 



; dx 



[l~2aa; + (a2-l)a;2]i/2- 



Now for a > 1, we take y — a — (a^ — l)a;, and 



/l (2/2-1)1/2 a2-l (a2_ 1)3/2 

For a < 1, we have y = a + {1 — a'^)x: 



a 



(1-2/2)1/2 l_^2 (1_ ^2)3/2 



arccos a. 



(Al) 



(A2) 



(A3) 



(A4) 



It is straightforward to verify that both expressions match smoothly at a = 1; /o(l) = V3- The next integral, 
Jo (a), can be reduced to the already calculated Iq. Taking x = , we get 



Ma) 



x^dx 



1 r {a^x^ + 2ax + 1 - 1 - 2ax)xdx 



(l + a2;)2(l-a;2)i/2 q,2 (1 + ax)2(l - j;2)1/2 



l-/n 



^dx 



2 a Jo (l + aa;)2(l -2;2)i/2 



1 + 3/o vr 2/o 



{a + x)dx 



Now it is useful to note that 



for all values of a. Thus, we arrive at 



/o (l + aa;)2(l-a;2)i/2- 
(a + x)dx _ 



Jo{a) 



(l + aa;)2(l-x2)i/2 

a2 - + 2 + (3a2 _ 2)/o(a) 



The remaining two integrals are easily calculated using the identity (p^), or 



Ii{a) — —a 



dlpja) 
da ' 



Ji (a) — —a 



c?Jo(a) 
da 



(A5) 
(A6) 

(A7) 

(A8) 



For convenience, all four integrals arc combined below 

1 



~ arccos a, [a < 1] 

I ^ 0.33333, [a = 1] 



(A9a) 
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Jo = < 



3^2 a(l+2a2) 

■(i_Q,2)2 + (i_a2)5/2 arccosa, 



^ ^ 0.26667, 



[a < 1] 
[a=l] 



' -a^+7ra2-|-2a-7r _ 3a2-2 

a3(l-a2) c«3(l-a2)3/2 



arccos a, 



f - TT ^ 0.19174, 



a3(l_a2 J 



[a < 1] 
[« = 1] 



(A9b) 



(A9c) 



f 2a^-SiTa^-na^+6TTa^+6a~3Tr , 12a'^-15a^+6 



Jl 



a3(l-a2)2 

9.6 - Svr ^ 0.17522, 



+ 



«3(l-a2)5/2 



arccos a, 



2a5_37i-a4_iiQ!3_)_g^Q,2_)_g^_3^ 12a4-l5Q.2-|-6 



a3(o2_i)2 



+ 



a3(a2-l) 



[a < 1] 

[a = l] 



^ln(a + (a2- 1)1/2), > 



(A9d) 



B. Orbit integration 

We describe here an accurate procedure to calculate orbits in a given potential. The four dimensionless 
potentials used in this paper and the corresponding angular momenta are given in Table ^ The angular 
momentum of an orbit of the given perigalactic distance and eccentricity is found by equating the orbital 
energy at the perigalacticon and the apogalacticon. 



We integrate the orbits using the following quadratures: 

dr 



r.2 / 2£-20(r) _ ^ 



1/2 



dr 



1/2 ' 



(Bla) 



(Bib) 



where r is the dimensionless radius, r = R/a or r = R/Rmax, and e is the dimensionless energy of the 
orbit. To avoid singularities at both the perigalacticon and the apogalacticon, we separate the integrals 
into two parts at some intermediate radius, rj_iy = (1 + e)rp. Then, we change the variables of integration 



individually for the two parts to remove siiiguls-rities. For Vp ^ r ^ '^div^ 

the new variable is ui = y a — 

for rdiv < r < Ta, it is U2 = \J ^ ^ "^TTe • Both parts of the integrals are regular in ui or U2 and can be 
calculated numerically with high accuracy. 



Note, that in the harmonic oscillator potential the orbits allow simple analytic solution. In Cartesian 
coordinates, 

X — Tp cost, y = Ta sini, (B2) 

where 

t = arctan i j-j-^ tan j . (B3) 



Table 1. 


Model 


Potential, ^(r) 


Angular momentum, j^(a, e) 


Hcrnquist model 
Harmonic oscillator 
Power-law, n = 1 


1 

l+r 

r= 3 
2 2 

r- 2 


a(l+a)(l+a + (l-Q)e) 
1— e 


Singular isothermal sphere 


Inr — 1 


-2 (1+e)^ 1 1+e 
" 2e 1-e 
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Fig. 1. — Integrals Iq, Jq, Ii, and Ji (eq. [A9|) for the Hernquist model. 
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Fig. 2. — Corrections to the velocity change relative to the point-mass approximation for the Hernquist 
model (cf eq. |^). The value of Xy is always close to zero, as it would be in a point-mass potential. For 
comparison, dashes show the fraction of the galactic mass inside the perigalacticon, xo = l^{Rp)- 
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Fig. 3. — The ratio of the correction factors, Xecc/Xst, for the eccentric and straight orbits of the satelhte. 
The hnes correspond to the orbits with eccentricities running from e = to 1 with a step of 0.1, in the 
Hernquist model for the galaxy. Here a = a/Rp, where a is the scale length of the model and Rp is the 
perigalacticon of the orbit. In the calculation of each curve, we assumed the same angular momentum of the 
straight orbit as for the corresponding eccentric orbit. Accordingly, the value of Xst may vary for each line. 
Adiabatic corrections are not included. Had they been, the correction factors for e — > would have been 
reduced dramatically. 



Fig. 4. — The correction factor Xecc for the power-law models (eq. [^) versus the eccentricity of the orbit 
(solid lines). For the isothermal sphere (n = 2), we also show the straight-path approximation, Xst (dots). 
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Fig. 5. — Time variation of the vertical component of the tidal force exerted on globular cluster NGC 6712 
by the Galactic bulge. Solid line is the actual force, whereas dashes show the Gaussian fit with the effective 
width T — 5.3tdyn- The force is in units of the N-body code. 
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Fig. 6. — The first and second order energy changes for globular cluster NGC 6712 (triangles). Both changes 
are plotted versus the initial energy of stars, Ei, and are normalized to the corresponding power of Ei. 
Error-bars indicate standard deviations in each energy bin. Here E{,ind is the initial binding energy of the 
system, and arrows are at the core and the half-mass radii. Solid lines are the analytic predictions (eq. [^). 
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Fig. 7. — Energy change of stars at the half-mass radius of NGC 6712 relative to the initial half-mass energy 
as a function of eccentricity of the orbit, e, and the corresponding perigalactic distance, Rp. The orbital 
energy is fixed at the observed value for the cluster (solid line) or at one- half and twice that value (dashed lines 
below and above the solid line, respectively). The star marks the eccentricity of NGC 6712 in our kinematic 
model (see §4.1). For illustration, we also show on the right-hand-side axis the approximate destruction time 
of the cluster due only to tidal shocks. The destruction time is defined as tdestr = Porb{Eh/ l^Eh), where 
Porb is the cluster orbital period, which varies with eccentricity. In contrast to Figure |^, which did not allow 
for adiabatic corrections, here tidal shocking vanishes as e ^ 0. 



